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Abstract 

An effective model in the isolated dimer limit of the Kitaev-type model on a decorated hon¬ 
eycomb lattice is investigated at finite temperature. The ground state of this model is exactly 
shown to be a chiral spin liquid with spontaneous breaking of time reversal symmetry. We elab¬ 
orate the finite-temperature phase diagram by using the mean-field approximation and Monte 
Carlo simulation. We find that the phase transition between the high-temperature paramag¬ 
netic phase and the low-temperature chiral spin liquid phase is always of second order in the 
Monte Carlo results, although a tricritical point appears in the mean-field phase diagram. The 
finite-size scaling analysis of the Monte Carlo data indicates that the phase transition belongs 
to the two-dimensional Ising universality class. 
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1 Introduction 

A Mott insulator without magnetic ordering down to the lowest temperature (T) has attracted 
considerable attention in strongly correlated electron systems [T]. This state, called the quantum 
spin liquid, is brought about by strong quantum fluctuations, which often become conspicuous in 
low-dimensional systems and in the presence of geometrical frustration of the lattice structures. 
Several candidates of quantum spin liquids have been found theoretically and experimentally in 
quantum spin systems on geometrically frustrated lattices. Among them, the chiral spin liquid 
(CSL), where the time reversal symmetry is broken despite the absence of magnetic order, is 
of particularly interesting and has been discussed, for instance, in the Heisenberg models on 
triangular and kagome lattices EHS]. 

As a simple realization of CSL, a Kiteav-type model defined on a decorated honeycomb 
lattice has been studied in the last decade mm- This is a variant of the Kitaev model on 
the honeycomb lattice [7|, for which the ground state is exactly shown to be a CSL. The 
decorated honeycomb lattice is constructed from the honeycomb lattice by replacing the sites 
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by triangles [see Fig. |TJa,)] . This model exhibits two types of CSLs, topologically trivial (with 
Abelian excitations) and nontrivial (with non-Abelian excitations), by changing the exchange 
constants [8] . The characteristics of the CSL ground states were also studied for the effective 
models in two limits, the isolated dimer limit and the isolated triangle limit [ 9 ]. Recently, the 
authors studied the finite-T properties of CSLs for the original Kitaev-type model by using 
the quantum Monte Carlo (MC) simulation [ 13 ]. In addition to the finite-T phase transition 
associated with the time reversal symmetry breaking, peculiar crossovers were found in the 
regions near the two limiting cases. Although the results were compared with the asymptotic 
behaviors expected in the two limits, the details were not reported in the previous study m- 
In order to further understand the mechanism of the crossovers and the phase transition, it is 
desired to elucidate the finite-T properties of the effective models in these limiting cases. 

In this paper, we investigate thermodynamic properties in the isolated dimer limit of the 
Kitaev-type model on the decorated honeycomb lattice. In this limit, the effective model is 
described by Ising-type Zi variables only. Performing the mean-field (MF) approximation and 
an unbiased MC simulation, we elaborate the finite-T phase diagram. In both calculations, we 
find that the system exhibits a phase transition associated with the time reversal symmetry 
breaking at finite T. The phase transition is always of second order in the MC simulation, 
whereas both first- and second-order transitions appear with a tricritical point in the MF 
approximation. We also find a crossover above the critical temperature, associated with the 
coherent alignment of local conserved quantities. Moreover, the analysis of the finite-size scaling 
indicates that the phase transition belongs to the two-dimensional Ising universality class. 

2 Model 

We begin with a variant of the Kitaev model defined on a decorated honeycomb lattice depicted 
in Fig. [D(a), whose Hamiltonian is given by [8] 

* = - £ E.V7W £ E-WW. (i) 

7 =x,y,z (ij ) 7 l=x,y,z (ij >; 

where aj is the y(= x,y,z) component of the Pauli matrices, representing an S = 1/2 spin 
located on site i; ( ij ) 1 and stand for the nearest neighbor (NN) pairs on 7 bonds within the 
triangles and on 7' bonds connecting the triangles, respectively [see Fig. [T[a)]. For simplicity, 
we take J x = J y = J z = J and J' x = J' = J' z = as in Ref. [8[. 

There are two kinds of conserved quantities in this model: W,i = II ied <J T an< ^ Wt = 
n, st cr,/’, where d ( t ) represents a set of the sites belonging to a dodecagon (triangle) loop on 
the decorated honeycomb lattice, and 7, corresponds to the bond component not included in 
the loop d or t among three NN bonds connected to the site i. These conserved quantities are 
Zi variables taking the values ± 1 . Note that W± changes its sign by the time reversal operation, 
as it consists of the product of three spins. The ground state of this model is a CSL in which 
the time reversal symmetry is broken by uniform alignment of Wt as well as Wd [8] . 

In the present study, we are particularly interested in the limit of J' J, where the dimers 
on the inter-triangle 7' bonds are weakly coupled by J. This is called the isolated dimer limit. 
When each dimer is regarded as a site, the interacting dimers form a kagome lattice, as shown in 
Fig. mb). The effective interaction between the dimers is derived by the perturbation expansion 
in terms of J/ J' [S]. The effective Hamiltonian is written in the form of 

H eS = T 0 - C h Y, W h - C t Y W h W tl W ta , (2) 

h (ht 1 t 2 ) 


2 


Phase Diagram of the Kitaev-type Model on . .. 


Nasu and Motome 




Figure 1: (a) Schematic picture of the model in Eq. JTJ on a decorated honeycomb lattice. 
Thin lines represent NN 7 bonds in triangles and thick lines represent NN 7 ' bonds connecting 
triangles, (b) A kagome lattice on which the effective model in Eq. (j2j is defined. 


where Wf, and W t are conserved quantities given by the projection of Wd and W t onto the 
dimer subspace; they are also Z 2 variables taking ± 1 . As shown in Fig. [ljb) , Wh and W t are 
defined on hexagons and triangles, respectively. Similar to W t , W t also changes its sign by 
the time reversal operation. The coefficients in Eq. © are given by Ch = ^ 77 f — rr 
and Ct = 2 §fg 7 77-; Eq is a constant, and the sum of (ht\t 2 ) is taken over three neighboring 
plaquettes [5]. 

In the following calculations, we parameterize the coefficients in Eq. © as Ch = x and 
C t = (1 — cc)/15 by introducing x £ [0 : 1]. The limit of J' J corresponds to x —> 1, where 
the second term in Eq. ([2]) is dominant over the third term. In the present study, however, 
we analyze the effective model in Eq. © in the whole region of x, i.e., 0 < x < 1, in order to 
elucidate the thermodynamic properties in this model. 


3 Method 

We investigate thermodynamic properties of the effective model in Eq. © using the MF ap¬ 
proximation and the MC simulation. In the MF approximation, the interaction between Wh 
and Wt in the last term in Eq. © is decoupled by introducing the mean fields (Wh) and (Wt). 
The resultant MF Hamiltonian per unit cell is given by 

nfi = -(x + (1 - x)(W t ) 2 )W h - 2(1 - x)(W h )(W t )W t + 2(1 - x)(W h )(W t ) 2 , (3) 

up to a constant. From this form, the self-consistent equations for the mean fields are given by 

(W h ) =tanh[/3{a; + (l-a ; )(Wt) 2 }], (4) 

(W t ) = tanh{2/3(l - x)(W h )(W t )}, (5) 

where /3 = 1/T is the inverse temperature (we set the Boltzmann constant fcs = 1). By solving 
the above equations and minimizing the free energy, we calculated the thermal averages (Wh) 
and (W t ). 

In the MC simulation, one can apply a classical MC method as the model in Eq. © includes 
only the Ising type variables Wh and Wt- We adopted the Wang-Landau algorithm [T4UT5] : 


3 





Phase Diagram of the Kitaev-type Model on . .. 


Nasu and Motome 



X 

Figure 2: Phase diagram of the effective model in Eq. ©• The horizontal axis x is defined so 
that Ch = x and Ct = (1 —x)/15 in Eq. (0. The single- and double-dashed green lines represent 
the critical temperature T c for the second- and first-order phase transition obtained by the MF 
calculation; the blue filled circle represents the tricritical point. The open circles connected by 
the solid red line and the open squares represent T c and T* obtained by the MC simulation, 
respectively. The dashed-dotted line indicates a crossover temperature T* estimated for an 
effective two-level system. See the text for details. 


after obtaining the density of states iteratively, we spent 10 6 MC steps for measurement and 
performed them 10 times, independently. We carried out the simulations for the finite-size 
clusters including N = 3 L 2 of Wh and Wt up to L = 50 under the periodic boundary conditions. 


4 Results 

4.1 Mean-field calculation 

The phase diagram obtained by the MF calculations is shown in Fig. [2] There is a phase 
transition between two phases, the high-T paramagnetic phase and the low-T CSL phase, 
associated with the time reversal symmetry breaking. The phase boundary is calculated as 
follows. As the CSL phase is characterized by the chiral order parameter (Wt), let us consider 
the expansion of the free energy per unit cell in terms of (Wt) as 

F MF = -i InTre-WZ* = f 0 + f 2 (W t } 2 + f 4 (W t ) 4 + ■■■ , (6) 

where the coefficients are given by 

/ 2 = 2(1 — :e){1 — 2/3(1 — x) tanh/3ir} tanh/te, (7) 

fi = 4/3(1 — a;) 2 [3{3 — 8/8(1 — x) tanh/3x}(l — tanh 2 /3x) + 8/3 2 (l — x) 2 tanh 4 f3x\. (8) 
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Figure 3: (a) T dependences of (a) the specific heat, (b) the entropy per site, (c) the mean 
square of the chiral order parameter k, and (d) the mean of Wh per hexagon at x = 0.0 obtained 
by the MC simulation, (e)-(f) Corresponding data at x = 0.8. For comparison, the specific 
heat for an effective two-level model, C v , is also shown in (e). The vertical dotted lines indicate 
T c and T*. The horizontal dotted lines in (b) and (f) represent S = - 11112 . 


When the transition is of second order, the critical temperature T c is determined by / 2 = 0 
with > 0, which leads to the equation 1 = 2/3(1 — a;) tanh/3a;. We find that this condition is 
indeed met in the large x region; the obtained T c is shown by the single-dashed green line in 
Fig. El In the vicinity of x = 1, T c arises linearly from zero, as the interaction coefficient Ct 
in Eq. (El) is proportional to (1 — x). While decreasing x, /a decreases to zero, which results 
in the tricritical point with / 2 = /4 = 0 at (x tr i, T tr i) ~ (0.615,0.596), as shown in Fig. El For 
x < a;tri, the phase transition becomes first order; in this region, T c is determined by comparing 
two local minima of F mf at (W t ) = 0 and {W t ) > 0. The obtained value of T c for the first-order 
phase transition is presented by the double-dashed green line in Fig. El 
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Figure 4: Scaling collapses for (k 2 ) at (a) x = 0.0 with T c = 0.4367, l/v = 1.027, and rj = 0.243 
at x = 0.0 and (b) x = 0.8 with T c = 0.1390, 1/v = 1.001, and 77 = 0.2545. 


4.2 Monte Carlo simulation 

We further examined the finite-T properties by the MC simulation, in which the effect of thermal 
fluctuations is taken into account beyond the MF approximation. Figure [5] show typical MC 
results for T dependences of physical quantities. First, we discuss the data at x = 0. FigureH a ) 
shows the specific heat for several cluster sizes. There is a peak that becomes sharper for larger 
L. The result indicates that the system exhibits a phase transition at T c ~ 0.44. Through this 
transition, all of the entropy In 2 is released, as shown in Fig.[3jb). Indeed, both W t and Wh grow 
coherently at T c . As shown in Fig. [3jc), the mean square of W t , ( k 2 ), where k = J2t Wt, 
abruptly increases below T c with decreasing T, indicating that the phase transition is associated 
with the time reversal symmetry breaking. At the same time, the mean of Wh per hexagon, 
(Wh) = j? J2h(Wh), also increases rapidly around T c , as shown in Fig. Hd). 

Next, we discuss the data at x = 0.8. In contrast to the case of x = 0, the specific heat 
shows two peaks, as shown in Fig. He). The low-T peak is similar to the one in Fig. H a )> 
which indicates a phase transition. On the other hand, the high-T peak is broad and does not 
show the system size dependence. This corresponds to a crossover, and we denote the peak 
temperature by T*. Figure Hf) shows the T dependence of the entropy per site. The entropy 
is released successively as decreasing T: one third of the entropy, |ln 2 , is released through 
the crossover at T* , and the rest two third | In 2 is released through the phase transition at 
T c . The latter is associated with W t , as indicated by the growth of (n 2 ) shown in Fig. H§)- 
On the other hand, the former comes from Wh', indeed, as shown in Fig. Hh), (Wh) increases 
around T* with decreasing T. Therefore, T* corresponds to the coherent growth of Wh, while 
T c corresponds to the time reversal symmetry breaking signaled by the abrupt growth of Wt- 
Thus, two conserved quantities develop at different T scales at x = 0.8, in contrast to the results 
at x = 0, where both grow simultaneously near T c . 

While the phase transition is driven by the last term in Eq. ©, the crossover at T* is 
caused by the second term. This is confirmed by considering the second term only, which gives 
an effective two-level system with the energy splitting 2x. The specific heat for this system 
is given by C v = /3 2 a; 2 / cosh 2 f3x. The result at x = 0.8 is presented by the black curve in 
Fig. He). This agrees with the MC result near and above T*. 
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To determine T c precisely and to clarify the nature of the phase transition, we performed the 
finite-size scaling for the order parameter (k 2 ). Figures Eta) andUJb) show the scaling collapses 
of (k 2 ) at x = 0 and x = 0.8, respectively. In both cases, all the data well collapse onto a single 
universal curve, which indicates that the transitions are of second order. Using the Bayesian 
scaling analysis for the optimization |16j . we obtain the critical temperature and the critical 
exponents as T c = 0.4367(1), 1/V = 1.027(9), and r] = 0.243(4) for x = 0 and T c = 0.13906(3), 
1/zz = 1.001(1), and r/ = 0.2545(8) for x = 0.8. The estimates of v and rj are close to those in 
the two-dimensional Ising universality class, v = 1 and r\ = 1/4. The scaling analysis indicates 
that the phase transitions belong to the two-dimensional Ising universality class. 

Carrying out the similar analyses while varying x, we obtain the MC phase diagram, as 
indicated by the red circles in Fig. [2] The x dependence of T c is similar to that in the MF 
result, while the values are reduced to roughly half due to thermal fluctuations. Interestingly, 
the scaling collapse works well for all x, indicating that the phase transition is always of second 
order in the MC results. This is in clear contrast to the MF result, where both first- and 
second-order transitions appear with the tricritical point. The results suggest that the tricritical 
behavior and the first-order transition are the artifact of the MF approximation. In addition 
to the phase transition, we show the crossover temperature T* by the orange squares in Fig. [2j 
which are estimated from the peak T in the specific heat. Near x = 1, T* well agree with 
T* ~ 0.8336x (dashed-dotted line) determined by the peak T of C v for the two-level system. 
With decreasing x from 1, the crossover at T* merges into the phase transition at T c , and Wh 
and W t grows simultaneously in the smaller x region. 

4.3 Comparison to the previous study 

Finally, let us compare the present results with those for the original Hamiltonian given in 
Eq. O, which were recently obtained by the authors ee The present results in the limit of 
x —> 1 well explain the asymptotic behaviors of T c and T* in the limit of J' J for the original 
model. Moreover, the scaling analysis for much larger system sizes compared to the previous 
study provides compelling evidence of the second-order phase transition belonging to the two- 
dimensional Ising universality class in the limit of J'/ J 3> 1. We note that tricritical behavior 
was found in the previous study EE and it looks qualitatively similar to the MF results in the 
present study. However, this will be just a coincidence, as the tricritical point in the original 
model is located far from the limit of J' J, close to the border between topologically-trivial 
and nontrivial CSL phases. This is far beyond the effective model used in the present analysis. 
Moreover, the tricritical behavior in the effective model is an artifact of the MF approximation 
as discussed above. Further study of more sophisticated effective models including higher-order 
perturbations will be necessary to clarify the origin of the tricritical behavior. This is left for 
future study. 


5 Conclusion 

In summary, we have investigated thermodynamic properties in the isolated dimer limit of the 
Kitaev-type model on the decorated honeycomb lattice. Using the MF approximation and 
the unbiased MC simulation, we found that the effective model exhibits a phase transition 
from paramagnet to the CSL at a finite T. The MC results indicate that the phase transition 
is always of second order, whereas the MF approximation gives both first- and second-order 
transitions with a tricritical point. The phase transition is associated with the time reversal 
symmetry breaking, which is identified by the conserved quantities Wt- In addition to the 
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transition, we found a crossover originating from the coherent growth of the other conserved 
quantities Wh- With increasing the inter-dimer coupling, the crossover is merged into the phase 
transition. Moreover, the finite-size scaling analysis shows that the phase transition belongs to 
the two-dimensional Ising universality class. We discussed the results in comparison with those 
for the original Kitaev-type model. 
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